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We present a directed unloading sand box type avalanche model, driven by slowly lowering the 
retaining wall at the bottom of the slope. The avalanche propagation in the two dimensional surface 
is related to the space-time configurations of one dimensional Kardar-Parisi-Zhang (KPZ) type 
interface growth dynamics. We express the scaling exponents for the avalanche cluster distributions 
into that framework. The numerical results agree closely with KPZ scaling, but not perfectly. 
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Avalanche phenomena are common in nature |jj . They 
are characterized by fast relaxation dynamics under a 
slow driving force. Models that describe such dynamics, 
e.g., so-called sandpile models, have been studied exten- 
sively for more than a decade following the work of Bak 
et al ||. Directed sandpile models are a special sub- 
class in which relaxation follows a directional rule ||, 
that is, the propagation of active sites occurs only in one 
direction and never backfires. The central issue in this 
type of research is whether the dynamics is critical, such 
that the avalanche distribution functions are scale invari- 
ant (power-law decay), and if so, whether these scaling 
properties are universal in the same sense as equilibrium 
critical phenomena. Our understanding of these issues 
is still restricted. There are only a few exactly soluble 
models, e.g., some deterministic Abelian directed sand- 
pile models ||], but most insight is still limited to nu- 
merical simulation data. The evidence for scaling and 
universality in other types of non-equilibrium dynamics 
is less ambiguous: In surface growth (another example 
of intrinsic critical behavior) several universality classes 
are well established; e.g., Edwards- Wilkinson (EW) 
and Kardar-Parisi-Zhang (KPZ) || type surface growth; 
population and catalysis type dynamics undergo absorb- 
ing state type phase transitions with distinct universality 
classes like directed percolation and directed Ising (|-|8). 

Efforts are under way to link avalanche dynamics 
to these better understood dynamic phenomena. This 
ranges from mappings to driven interfaces Jsj — and di- 
rected percolation p"2| , |l3| ] to using the concepts of renor- 
malization [§~4|,[l5| and universality classes [|6jl7|l^]. In 
this letter we introduce a two dimensional (2D) directed 
avalanche model in which the scaling exponents of the 
avalanche dynamics are directly related to the exponents 
of KPZ growth of a ID interface. 

The physical system we have in mind is a sand box 
with a movable retaining wall to let out sand from the 
bottom of the slope, see Fig. ([!]). The retaining wall is 
lowered very slowly, such that grains tumble out sporad- 
ically forming distinct avalanches instead of a continuous 
flow. We model the sand surface by continuous height 
variables h(x,y), with respect to a 2D square lattice, 



which is rotated over 45°, meaning that in the even (odd) 
y rows x takes only even (odd) integer values. The re- 
taining wall is placed at the y — boundary. The slope is 
stabilized by the following constraint. The surface par- 
ticle in column [x, y) is supported by the two columns 
(x ± 1, y — 1) just below it, and must be lower than the 
lowest of the two increased by an amount s c , 

h(x, y) < min [ft (a; + 1, y - 1), h(x - 1, y - 1)] + s c . (1) 

An avalanche is triggered by selecting the highest site, 
(xi,0), at the y = wall boundary (it is the i-th 
avalanche) and reducing its height by a random amount, 
< rji < s c . This represents the lowering of the retaining 
wall. Next, all sites that violate the stability condition 
topple according to the rule, 



h(x, y) — > min [h(x + 1, y — 1), h{x 

+ Vi(x,y), 



i,y-i)] 



(2) 



where < rji(x,y) < s c are uncorrelated random num- 
bers. This toppling continues until the whole system is 
stable again. Since the toppling of a site in row y can 
only effect the stability of two sites immediately above it 
in row y + 1, the sites can be updated row by row starting 
from the y = boundary. 

This process is idealized compared to a real unloading 
sand box in the sense that the toppled grains drop out 
without disturbing the already stabilized lower regions of 
the surface. It is possible to justify this as the low gravity 
or strong bond limit where the falling sand does not gain 
enough momentum to disturb the stabilized surface on 
its way out. Without this idealization we loose the di- 
rected nature of the dynamics and the following surface 
growth interpretation. 

The row-by-row toppling sequence (j^) can be reinter- 
preted as a dynamic rule for a ID growing interface, in 
which the y coordinate plays the role of time. Each 
stable surface configuration represents a space-time con- 
figuration of the ID interface. Imagine creating an ini- 
tial stable surface configuration, before the retaining wall 
starts to drop: choose an arbitrary configuration with all 
< h(x, 0) < s c in row y = 0. Next, apply Eq.(||) to 
all sites in the next row, y = 1, to create the next slice 
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of the surface. Repeat this for all higher rows. The con- 
figuration in every row is like a ID interface evolving in 
time t = y. Fig.|] illustrates how this interface prop- 
agates during each time step, y — > y + 1. The upper 
panel shows the first half of the update (the drawn to 
the dashed line). This is the deterministic part of the 
propagation (the min[] operator) in Eq.(|^). The lower 
panel illustrates the second half of the update, where the 
heights increase by a random amount < rj < s c . The 
first step removes material, and the second step deposits 
particles. 

This type of interface dynamics almost certainly be- 
longs to the KPZ universality class. Eq.(||) can be rewrit- 
ten as 

h(x, t) = - [h(x + l,t-l) + h(x -l,t- 1)] 



- - \h[x + l,t - 1) - h(x - l,t - 1) 



(3) 



which is a discrete version of the KPZ Langevin equa- 
tion @, 



hi 



To be absolutely sure, and also to make sure that A is 
large enough such that corrections to scaling from the 
EW point (at A = 0) do not interfere, we checked numer- 
ically the behavior of the surface width A(L,t), defined 

as A 2 = (^(h - {h)) 2 y The scaling properties of KPZ 
growth are known exactly in 1+1 dimensions. Starting 
from a flat surface at y = 0, the width increases as 
for < t < L z and saturates at A ~ L a for t > L z , 
with L the lattice size in the x direction. The exponents 
are exactly equal to a = ^, (3 = | and z = a/ [3 = |. Our 
numerical results are shown in Fig.([|) and are consistent 
with the KPZ values. 

Throughout this paper we present numerical results in 
terms of finite size scaling plots of effective exponents, 
like a(L) in Fig.([|)(a). Global straight line fits to log- log 
plots of, e.g., A ~ L a are inaccurate. a(L) is obtained 
from fitting the form A ~ L a to two nearby values of L. 
The approach to L — > oo in Fig.(^) is consistent with a 
leading corrections to scaling exponent y = — |. 

The characteristic feature of self-organized criticality 
(SOC) is the lack of typical avalanche length, width, 
depth, or mass scales. The probability distributions fol- 
low power-laws, like P w ~ w~ Tw characterized by the 
scaling exponents, r;, t w , Tg and r m . Our numerical 
simulation results confirm the existence of scale invari- 
ance. The critical exponents converge well, see Fig.(^). 
The length I is the maximum y coordinate the avalanche 
reaches. The width w is the maximum departure in the x 
direction, \x — Xi\, from the trigger point n. The depth 
8 is the maximum height change, hi — caused by 



the avalanche, and the mass is the total amount of sand 
carried off by the avalanche. 

The meta-distribution function, P(l, w, S), should obey 
the scaling form 



P(l,w,S) = b- a P(b- z l,b- 1 w,b- a 6). 



(•5) 



with b an arbitrary scale parameter. The exponents, a, 
z and a, are expected to be robust with respect to de- 
tails of the dynamic rule, and thus define the universality 
class. Single parameter distributions, like P w ~ w~ Tw , 
follow by integrating out the other two parameters. This 
implies the following expressions for the r exponents 



cr — 1 — a o — 1 — z 
n = , t w = a - z - a, Tg = . 



or inverted, 



T w ~l 



T w -l 
TS — 1 



a = t w + z + a. 



(6) 



(7) 



In Fig. (g) we replot the numerical finite size scaling esti- 
mates of the r exponents in terms of a, z, and a. 

The values z = 1.52±0.02 and a = 0.46±0.01 are very 
close to those of ID KPZ growth. Every stable slope 
configuration represents a possible world sheet of a ID 
KPZ type interface, and every avalanche the difference 
between two such world sheets. Therefore it is natural 
to expect that the length (depth) of the avalanche scales 
with the KPZ value for z (a). 

The distribution of avalanche cluster sizes is the most 
commonly studied and experimentally the most accessi- 
ble property of SOC. Its exponent is linked to a and z in 
the following manner. At the start of the avalanche, the 
height of a boundary site (y = 0) is lowered on average 
by \s c . Thus, for a sandbox of width, L Xl the boundary 
row is lowered by ^s c after L x avalanches. In the sta- 
tionary state, the entire surface moves down on average 
by \s c and the average amount of removed sand is equal 
to L y x ^L x s c . Thus, the average mass of an avalanche 
must be equal to 



(to) = J dmmP m (m) — —s c L y . 



(8) 



This is analogous to conservation of current in conven- 
tional deposition type avalanche systems (see, e.g., 
Assume that the avalanche clusters are compact, i.e., that 
the sizes of holes of unaffected regions inside an avalanche 
do not scale with the avalanche size. In that case, the 
mass scales as m ~ / x w x 6, and we can use the meta- 
distribution function to evaluate Eq.(||), as 

pL y p oo p oo 

(to) ~ / dl dw dS lw8 P(l,w,S) 
Jo Jo Jo 

/• oo /• oo /> oo 

+ m Ly / dl dw d8 P(l,w,8) (9) 
Jl v Jo Jo 
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This applies when the box is wide and deep enough that 
L y is the only limiting factor to the avalanches. The 
first term accounts for all avalanches that fit inside the 
box, and the second term for the ones that reach the L y 
edge and thus are prematurely terminated. The first in- 
tegral scales as L y a + 2 + 2z + 2a )/ z f or j ar g e while the 



in, e.g., Ref. j 
effect like this. 



L9| is not sufficient to detect a small 



second one scales as L 
m ~ I x w x S so rnr 



(-cr+l+z+a)/z 



-(l+z+a)/z 



U 



Eq.(g) scale in the same way, as 
Solving Eqs.(|) and © for a gives 



We have assumed 
The two terms in 

(10) 



a = 2 



2a 



(11) 



Numerically we find a — 4.43±0.05, see Fig.(||), in agree- 
ment with the numerical values for a and z, and also with 
their KPZ values (which yields a — 9/2). 

In conclusion, we introduced a directed unloading sand 
box model, in which the stable slope configurations obey 
KPZ type scaling, and the avalanches represent the differ- 
ence between two such KPZ world sheets. The avalanche 
distribution exponents can be reformulated in KPZ lan- 
guage, Eq. (]?]), and the numerical results agree with the 
KPZ values. 

An intriguing issue for further study is whether the 
slight systematic deviations from the KPZ values in 
Fig.([|) are for real or just an artifact of, e.g., finite size 
scaling (the avalanches cover only a small part of the 
surface). In KPZ dynamics the scaling exponents follow 
from an ensemble average over independent MC runs, 
i.e., over a large set of independent world sheets. The 
avalanche dynamics performs this ensemble average in a 
correlated manner. All subsequent world sheets are iden- 
tical except for the avalanche area. Correlated MC runs 
like this might well change the scaling exponents (but 
only slightly apparently) , and this might prove a key fea- 
ture of avalanche type SOC. The dashed curve in Fig.|](b) 
shows initial results for the global surface roughness in 
such avalanche correlated MC runs. The surface is no- 
ticeably rougher in amplitude but the finite size scaling 
values of the exponent /3 are somewhat smaller (consis- 
tent with the values for a and z in Fig.||). Whether they 
will converge in the end to a value (3 < \ remains to be 
seen. 

In this context it is noteworthy that Eq.(|^) applies 
also to recent solutions to the stochastic directed sand- 
pile model (SDSM) of Paczuski et.al. Jisj and Kloster 
et.al. [|l9| where the avalanche dynamics is related to 
Edwards- Wilkinson (EW) M surface growth, with z = 2 
and a = 1/2; and also that the exactly soluble directed 
sandpile model (DSM) of Dhar and Ramaswamy satisfies 
Eq.(^) with z — 2 and a — 0. However, that does not 
resolve the issue, because both models are intrinsically 
simpler than ours, and the numerical accuracy of SDSM 



This research is supported by the National Science 
Foundation under grant DMR-9985806. 



[1 
[2 
[3 
[4 
[5 
[6 

[-: 

[8 

[9 
[10 

[11 

[12 
[13 

[14 

[IS 

[16 

[17; 

[is: 

[19 

[20 



P. Bak, How nature works: the science of self-organized 
criticality (Copernicus, New York, NY, USA, 1996). 
P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. Lett. 
59, 381 (1987). 

D. Dhar and R. Ramaswamy, Phys. Rev. Lett. 63, 1659 
(1989). 

S. F. Edwards and D. R. Wilkinson, Proc. R. Soc. Lon- 
don A 381, 17 (1982). 

M. Kardar, G. Parisi, and Y.-C. Zhang, Phys. Rev. Lett. 
56, 889 (1986). 

W. Kinzel, Ann. Isr. Phys. Soc. 5, 425 (1983). 
For a review see J. Marro and R. Dickman, Nonequi- 
librium Phase Transitions in Lattice Models (Cambridge 
University Press, Cambridge, 1996). 

J. D. Noh, H. Park, and M. den Nijs, Phys. Rev. E 59, 
194 (1999). 

K. Sneppen, Phys. Rev. Lett. 69, 3539 (1992). 

M. Paczuski and S. Boettcher, Phys. Rev. Lett. 77, 111 

(1996). 

K. B. Lauritsen and M. J. Alava, arXiv.org e-print 
(1999), [sond-mat/9903346| . 

B. Tadic and D. Dhar, Phys. Rev. Lett. 79, 1519 (1997). 
A. Vazquez and O. Sotolongo Costa, J. Phys. A 32, 2633 

(1999) . 

L. Pietronero, A. Vespignani, and S. Zapperi, Phys. Rev. 
Lett. 72, 1690 (1994). 

A. Vespignani, S. Zapperi, and L. Pietronero, Phys. Rev. 
E 51, 1711 (1995). 

A. Ben-Hur and O. Biham, Phys. Rev. E 53, R1317 
(1996). 

A. Mehta, J. M. Luck, and R. J. Needs, Phys. Rev. E 53, 
92 (1996). 

M. Paczuski and K. E. Bassler, Phys. Rev. E 62, 5347 

(2000) . 

M. Kloster, S. Maslov, and C. Tang, Phys. Rev. E 63, 
026111 (2001). 

R. Pastor-Satorras and A. Vespignani, Phys. Rev. E 62, 
6195 (2000). 



3 




X 1/w 1/w 




FIG. 2. KPZ type growth dynamics 




FIG. 3. Monte Carlo (MC) results for the global interface 
width: (a) finite size (L x ) estimates for the saturated surface 
width exponent a; (b) finite time estimates for the transient 
surface width width exponent /3 from a flat initial configura- 
tion. The drawn line represent conventional MC simulations 
with an uncorrelated world sheet ensemble, and the dashed 
line are data from avalanche correlated MC runs. 4 



